Nonequilbrium-induced metal-superconductor quantum phase transition in graphene 
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We study the effects of dissipation and time-independent nonequilibrium drive on an open super- 
conducting graphene. In particular, we investigate how dissipation and nonequilibrium effects mod- 
ify the semi-metal-BCS quantum phase transition that occurs at half-filling in equilibrium graphene 
with attractive interactions. Our system consists of a graphene sheet sandwiched by two semi-infinite 
three-dimensional Fermi liquid reservoirs, which act both as a particle pump/sink and a source of 
decoherence. A steady-state charge current is established in the system by equilibrating the two 
reservoirs at different, but constant, chemical potentials. The graphene sheet is described using 
the attractive Hubbard model in which the interaction is decoupled in the s-wave channel. The 
nonequilibrium BCS superconductivity in graphene is formulated using the Keldysh path integral 
formalism, and we obtain generalized gap and number density equations valid for both zero and finite 
voltages. The behaviour of the gap is discussed as a function of both attractive interaction strength 
and electron densities for various graphene-reservoir couplings and voltages. We discuss how trac- 
ing out the dissipative environment (with or without voltage) leads to decoherence of Cooper pairs 
in the graphene sheet, hence to a general suppression of the gap order parameter at all densities. 
For weak enough attractive interactions we show that the gap vanishes even for electron densities 
away from half-filling, and illustrate the possibility of a dissipation-induced metal-superconductor 
quantum phase transition. We find that the application of small voltages does not alter the essential 
features of the gap as compared to the case when the system is subject to dissipation alone (i.e. zero 
voltage). The possibility of tuning the system through the metal-superconductor quantum critical 
point using voltage is presented. 

PACS numbers: 03.65.Yz, 64.70.Tg, 74.78.-w 



I. INTRODUCTION 

The landmark experimental realization of an isolated 
graphite monolayer, or graphenei^, has sparked intense 
theoretical and experimental interest in the material over 
the last few years 3 -. A source of interest in the study 
of graphene is the unique properties of its charge car- 
riers. At low energies, these charge carriers mimic rel- 
ativistic particles, and are most naturally described by 
the (2+l)-dimensional Dirac equation with an effective 
speed of light, c ~ vf ~ 10 6 ms _1 . The fact that 
graphene is an excellent condensed-matter analogue of 
(2+l)-dimensional quantum electrodynamics (QED) has 
been known to theorists for over 20 yearsi££. However, it 
was not until the spectacular experimental realization of 
isolated graphene that experimentalists began observing 
signatures of the QED-like spectrum in their laboratories. 
Consequences of graphene's unique electronic properties 
have been revealed in the context of anomalous integer 
quantum-Hall effect^ and minimum quantum conduc- 
tivity in the limit of vanishing carrier concentrations^. 

In addition to its importance in fundamental physics, 
graphene is expected to make a significant impact in 
the world of nano-scale electronics. Research efforts 
in developing graphene-based electronics have been fu- 
eled by a strong anticipation that it may supplement 
the silicon-based technology which is nearing its limits 3 . 
Graphene is a promising material for future nanoelec- 
tronics because of its exceptional carrier mobility which 
remains robustly high for a large range of temperatures, 
electric-field-induced concentrationsii 2 ^^, and chemical 



doping^. Indeed, recent experiments have explored the 
possibilities of in-plane graphene hctcrostructures by en- 
gineering arbitrary spatial density variation using local 
gate o 1 ^ 11 ' 12 . The application of local gate techniques to 
graphene marks an important first step on the road to- 
wards graphene-based electronics. 

From a theoretical point of view realizing graphene 
nanoelectronics requires a theoretical understanding of 
open nonequilibrium graphene. Naturally, graphene in 
nano-circuits is subject to decoherence effects due to its 
coupling to external leads via tunnel junctions. Further- 
more, a nonequilibrium treatment of graphene becomes 
necessary when a charge current is driven through it. 
To this date, effects of dissipation and nonequilibrium 
drive on graphene electronic properties have not been 
addressed. The focus of this paper is to show a theoreti- 
cal framework in which these effects can be studied and 
illustrate how they give rise to striking influences on the 
equilibrium properties of graphene. 

This work considers dissipation and nonequilibrium 
effects on superconducting graphene. Besides the pos- 
sibility of superconductivity in graphene by proxim- 
ity effect^, some works suggested the potential of 
achieving plasmon-mediated singlet superconductivity in 
graphen o 14 ' 15 . Several groups have investigated the 
equilibrium mean-field theory of superconductivity in 
graphene using the attractive Hubbard model on the 
honeycomb lattice. Uchoa and Castro-Netc 14 studied 
spin singlet superconductivity in graphene at various fill- 
ings by considering both the usual s-wave pairing as 
well as pairing with p + ip orbital symmetry permit- 
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ted by graphene's honeycomb lattice structure. Zhao 
and Paramekantii^ examined the possibility of s-wave 
superconductivity on the honeycomb lattice. Both works 
show that (in the absence of p-wave pairing) half-filled 
graphene displays a semimetal-superconductor quantum 
critical point at a finite critical attractive interaction 
strength u c . Away from half-filling, the system exhibits 
Cooper instability at any finite u and thus undergoes the 
usual BCS-BEC crossover as u is increased. The difficulty 
in achieving superconductivity at half-filling is a result of 
the vanishing density of states at the Dirac point and the 
absence of electron screening. 

In this work, the superconducting graphene sheet is 
subjected to dissipation and noncquilibrium drive by cou- 
pling it to two semi-infinite particle reservoirs via tun- 
nel junctions. The geometry of the system is shown in 
FigQ] While the two reservoirs are independently held 
in thermal and chemical equilibrium at all times, an out- 
of-plane steady-state current through graphene is estab- 
lished by equilibrating the reservoirs at two different, 
but constant, chemical potentials. The leads act as infi- 
nite reservoirs and are assumed to be held at a common 
temperature T at all times. Nonequilibrium theory of 
BCS superconductivity is formulated using the Keldysh 
path integral formalism, and the resulting nonequilib- 
rium mean-field equations are used to investigate the gap 
behaviour at and near half-filling for various attractive 
interaction strengths. The gap is plotted in the parame- 
ter space of filling n and the interaction strength u (see 
FigjS]), and our results can be directly compared to the 
gap phase diagram in Fig. 2 of the work by Zhao and 
Paramekantii£. 

Our main results are now qualitatively summarized. 
We find that the gap is generally suppressed in the pres- 
ence of leads. As the paper will discuss in detail, the key 
to understanding our findings is to notice that the dissi- 
pation of electrons into the leads act as a pair-breaking 
mechanism for the Cooper pairs in the central graphene 
sheet. This mechanism, and hence the suppression, is 
present at both zero and finite voltages and for all elec- 
tron densities. As a consequence, the Fermi liquid ground 
state of the system remains stable against Cooper pairing 
up to some density-dependent finite attractive interac- 
tion strength u c (n) at all densities. With respect to the 
gap phase diagram, dissipation gives rise to a finite region 
around half-filling in which the gap vanishes (see FigjS]). 
From these results, we infer that dissipation induces a 
metal-superconductor quantum phase transition at all 
fillings, for which the tuning parameter is the attractive 
interaction strength u. The qualitative behaviour of the 
gap is not appreciably different in the zero and finite 
voltage cases as long as the voltage is small, i.e. V <C T, 
where T denotes the average tunneling rate of electrons 
between graphene and the two leads. Finite voltage mod- 
ifications, however, result due to voltage-induced changes 
in the graphene electron density. 

The paper is organized as follows. In SecHH we in- 
troduce the Hamiltonian which models our heterostruc- 
ture. The mean-field treatment of the model is formu- 



lated on the Keldysh contour in Sec lIIIl In Sec lIII Bl the 
nonequilibrium gap and number density equations will 
be derived. The results are presented in Sec lIVI The ef- 
fects of dissipation in the absence of voltage is discussed 
in Sec lIV Al while the finite voltage effects are included 
in Sec lTVBl We conclude in Sec|V] 

II. MODEL 

The lead-graphene-lead heterostructure considered in 
this work is shown in FigQ] Graphene is located on the 
z = plane, and each of its sites is labeled using two 
coordinates r ; = (xi , yi , Zi = 0) . The semi- infinite metal- 
lic leads extend from both sides of the graphene sheet 
for z > and z < 0. We assume the leads are sep- 
arated from graphene by thin insulating barriers, and 
the tunneling of electrons through each of the barriers 
can be described by phenomenological tunneling param- 
eters. Full translational symmetry is present along the 
planes parallel to the xy-plane for z ^ while only the 
discrete translational symmetry of the graphene lattice 
is present at z — 0. The leads are assumed to be in 
thermal equilibrium with their continuum of states occu- 
pied according to the Fermi-Dirac distribution, f a {^>) 
I + cxp ([3(uj — Ha))]~ , where a = L(left), i?(right) la- 
bels the leads. An electric potential bias is set up in the 
out-of-plane direction by tuning the chemical potentials 
of the leads to different values. 

left lead /i£ Honeycomb lattice layer 




FIG. 1: A schematic of the system considered. Chemical po- 
tential mismatch in the two leads will lead to a charge current 
parallel to the z-axis. 



The Hamiltonian consists of three parts, 

H Hsys Hres -(- H sy S — reS . (1) 

The central graphene sheet is modeled using the attrac- 
tive Hubbard model on the honeycomb lattice. The 
kinetic term is a tight-binding description for the tt- 
orbitals of carbon that includes nearest- and next- 
nearest- neighbour hopping processes. The on-site inter- 
action strength is parametrised by U . The Hamiltonian 
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for the layer is 

H sys = -t ( c l* c j,v + h - c -) 

-t' ( C i,<r C J> + h - c -) ~ U J2 C h C ll C *'t Ci ^- ( 2 ) 

cj a (cj.cr) creates (annihilates) electrons on site rt of the 
graphene honeycomb lattice with spin a (a =T,i)- U 
is assumed positive due to attractive interaction, and t 
and t' are the nearest- and next-nearest neighbour hop- 
ping parameters, respectively. Specific values for t and t' 
have been estimated^ by comparing a tight-binding de- 
scription to first-principle calculations. Following their 
estimates, we take t = 2.7eV and fix t' ft = 0.04. 

The honeycomb lattice can be described in terms of 
two inter-penetrating triangular sublattices, A and B 
(see Fig|2]). Each unit cell is composed of two atoms, 
one each of type A and type B. Primitive translation 
vectors, e! and e 2 , are 

ei = (V3,0) e 2 = (- V3/2, 3/2) e 3 = ei +e2, (3) 

where they are expressed in units of a, the distance be- 
tween two nearest carbon atoms. Any A atom is con- 
nected to its nearest neighbours on the B lattice by three 
vectors 



di = (0,1) 

d 2 = (-x/3/2,-1/2) 

d 3 = (^3/2, -1/2). 



(4) 




FIG. 2: Graphene honeycomb lattice, ei and are the unit- 
cell basis vectors of graphene with lattice constant \^3a w 
2.46 A (a ss 1.42^4). A unit cell contains two carbon atoms be- 
longing to the two sublattices A (white circles) and B (black 
circles) . All nearest- and next-nearest-neighbour hopping ma- 
trix elements are —t and —t', respectively. 
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FIG. 3: A diagram illustrating the type of tunneling process 
that are considered in this work. The diagram is an edge-on 
view of the interface between the graphene sheet and a lead. 
The only tunneling events that are allowed are those in which 
the (x,y) coordinates of electrons remain unaltered. Thus, 
while the lower two processes in the diagram are allowed, 
tunneling of the type shown at the top is disallowed. 



In momentum space, the kinetic term reads 



I/ K _ 

sys 



k,CT 



A k 3k 

.9k A k 



ak,, 



where 



3k = 



-t' e * k ' ei + c - c 



(5) 

(6) 
(7) 



Components of the pseudospinor, a£. a and a , describe 
quasiparticles that belong to sublattice A and B, respec- 
tively. Here, denotes the number of lattice sites in 
a triangular sublattice. N = 2N& will denote the total 
number of sites on the honeycomb lattice. 

Coupling between leads and the graphene sheet is mod- 
eled using the following Hamiltonian, 



J T D A ~ 



a=L,R i,a 



C, a is a phcnomcnological tunneling matrix that de- 
scribes the tunneling of an electron between site i on 
the graphene sheet and an adjacent site on lead a (see 
FigOJ) ■ We only consider lead-graphene tunneling pro- 
cesses in which (x, y) coordinates of the electron in the 
initial and final states are the same. This assumption 
simplifies various computational steps without altering 
the qualitative features of the final results. C\ a ak cre- 
ates an electron in lead a at coordinates (xi, yi) with spin 
a and longitudinal momentum k z . We assume here that 
the tunneling parameters are independent of frequency 
and momentum but maintain their lead dependence in 
order to describe possible asymmetries in the lead-layer 
couplings. In momentum space, the tunneling Hamilto- 
nian EqJS] becomes 



H, 



sys— res 



E, f dk z 1 s 

a k.cr 



A 



k,k z 



5k,<7 



,). (9) 



The inplane momentum, k, is the component of momen- 
tum parallel to the graphene plane and the out-of-plane 



(8) momentum, k z , is its component normal to the plane. 

) corresponds to an electron mode 



At 



k,fc 2 
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propagating in "sublattice A(B)" in lead a with spin 
a and wavevector k. Although the full inplane trans- 
lational symmetry of the leads implies that k can take 
on any value in R 2 the tunneling assumption (see Fig[3]) 
tells us that the only modes that tunnel are those with k 
values that are the allowed modes of the triangular sub- 
lattices in the graphene sheet. All other inconsequential 
modes can eventually be integrated out in the path in- 
tegral sense and will merely contribute a multiplicative 
factor in front of the partition function. Therefore, we 
will not consider these modes further. 

Both leads are assumed to be Fermi liquids 

O.A J k,fc a ,<7 

£k,fc, ( A i,kz,a,a A Kk z ,<T,a + B l,k^a, a B k,k z ,a,c^j , (10) 



with a separable dispersion 



I 2 kl 



£k,fc z = Ck + £fc z = 



2m e 2m e 



(11) 



Besides their role as a particle pump/sink, the leads play 
an important role as a heat sink. An important assump- 
tion we make is that any heat generated in the interacting 
region due to the application of a transverse electric field 
is efficiently dissipated into the leads so as to prevent 
build up of heat in the region. This is a well-justified 
assumption because the leads are assumed to be infinite 
and the interacting region has a thin profile. 

In equilibrium (fj, res — (Ar — /ir,), the central system is 
expected to reach chemical equilibrium with the reser- 
voirs in the long-time limit so that fi sys — fi res - In 
the out-of-equilibrium case, the system is coupled to two 
reservoirs that are not in chemical equilibrium. There- 
fore, although the electron distribution in the interacting 
system reaches a static form in the long-time limit, it 
is in no way expected to have an equilibrium form due 
to constant influx (outflux) of particles from (into) the 
leads. 



III. KELDYSH PATH INTEGRAL 
FORMULATION 

In this section, we formulate a theory of nonequi- 
librium BCS superconductivity in graphene using the 
Keldysh functional integral formalism. The theory is first 
expressed in terms of a Keldysh partition function using 
coherent states of fields defined on the time-loop Keldysh 
contour, C. Following a Hubbard-Stratonovic decoupling 
of the quartic interaction term in the pair channel a BCS 
theory for superconducting graphene is obtained by as- 
suming a static, homogeneous gap, integrating out both 
leads and graphene electrons, and extremizing the effec- 
tive action with respect to the gap. The resulting mean- 
field equations, which are a noncquilibrium generaliza- 
tion of the corresponding equilibrium equations^, are 
analyzed in the remainder of the paper. 



The starting Keldysh generating functional reads 

Z K = J 3>{a, a, b, b, A, A, B, B}e isK , (12) 

where 

(13) 



sys ' res ' sys — res ' 

If we introduce 4-component spinors defined in Nambu- 
sublattice space for both graphene electrons and leads 
electrons 



$k,k*,c«(i) 



a k,t(£) 

= [ a_ k ,|W 
&m(*) 

5-k,i(*) 

-<4k,k z ,T,a(t) 
— ^-k,-k,,J.,a(*) 

-£>k,k z ,T,a(i) 
-■B-k,-k»,i,a(*) 



(14) 



(15) 



the actions in Eq[T51 become 

x [idt - A k rf - 3k rf - &r? 4] ^(t) 
+ U dty2[a^(t)a td (t)a l . l (t)a hl (t) 

+M*)5i,^*)M*)M*)] . ( 16 ) 



X [IC 



d t - Ck,fc.7f) $k,k.,a(t), ( 17 ) 



and 



S. 



K 



dt 



dk. 



2n £ Ca iv A ? 



[*k,I^,a(*)^k(*) + 0k(*)*k,k.,a(t)] • (18) 

t± are 2x2 matrices given by 



4 = -«±ir;), 



(19) 



where t% y z are Pauli matrices. Superscript v indicates 
the space in which the matrices act; A (N) denotes sub- 
lattice (Nambu) space. The quartic interaction term 
in EqfTH] is decoupled using Hubbard-Stratonovic fields 
Af(t) and Af (i). In the BCS mean- field approxima- 
tion, where this field is assumed static and homogeneous 
(i.e. Af(t) = Af (t) = A), the resulting action of the 
system reads 



[•a \ Af NA * Af A 

+UAt% + C/AV5] Mt) - 2C/ 1 A| 2 . (20) 
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The self-consistency condition for the gap is 

A = {a U a iA ) (t) = {b lA h A ) (t). (21) 

The time-loop contour integral is carried out by first 
splitting every field into two components, labeled "+" 
and "— ", which reside on the forward and the backward 
parts of the time contour, respectivel y 18 ! 19 ' 20 . The con- 
tinuous action then becomes 



S K = 



dt[jSf+(t)-jS?_(f)], 



(22) 



where „Sf± (t) is the Lagragian corresponding to the action 
defined in Eq[l3] written in terms of + (— ) fields. When 
time-ordered products of Heisenberg fields in the theory 
are constructed on the Keldysh contour we obtain four 
Green functions 



iG T (t,t') 
iG f (t,t') 
iG < (t,t') 
iG > (t,t') 



<T + (t)? + (f)) 
(T_(t)T_(f)> 
<T + (t)T_(f)) 
(T_(i)T + (t')> 



Because these Green functions are not linearly inde- 
pendent, a linear transformation of the fields from the 
Kadanoff-Baym basis (+,— ) to the Keldysh basis (cl,q 
for bosons; 1,2 for fermions) is commonly performed. For 
bosons, the barred fields are related to the unbarred fields 
simply by complex conjugation and thus the transforma- 
tion is identical for both, 



T„ 



1 

V2 



1 1 

1 -1 



(23) 



For fermions, unbarred fields transformed in the same 
manner as Eq[23] For barred fields, we choose a different 
transformation 18 



Ti 
T 2 



J_/l -1 

vn 1 1 



(24) 



In order to express the Keldysh action, Eqf221 in 
the Keldysh basis it is now appropriate to define 8- 
component spinors for graphene electrons and leads elec- 
trons defined in the Nambu-sublattice-Keldysh space. 
Since we are interested in steady-state properties of the 
system, it is useful to first Fourier transform the fields 
into frequency space. We define the 8-component spinors 
as 



( a M 


\ 


a -k,l 




b h 




b ii 




a h 




a -k,i 




b h 






/ 



/ 


A 1 

^k,k z ,T,a 


\ 




A 1 






pi 

_ k,k z ,\,a 






~B-k,-k z ,l,a 






A 2 






A 2 






R 2 

D k,k z A,u 




V 


R 2 

D -k.-k z ,l,a 


1 



(25) 



where k = (k, ui) is the energy-momentum 3-vector. The 
action (Eq [2^|) then becomes 



3 K 

J sys 



+ U [A qT ? + A*t. 



9^rt 



A c/ r 



(A ci < + A c >-)rf]}^ 
- 2U [A*,A g + A*A C /] , (26) 



^res 



dkz 
2tt 



J2^k,k z , a 



x {(^(fc)rf 
(^(fc)rf - 



N T K 



(27) 



and 



sys— res 



-^r/Z^ [**,*„aV'fc + ^k,k z , a ] ■ (28) 



2tt 



Here, J k = j^—YlkJ ^7' an d r T,l are 2x2 matrices de- 
fined by 



T T,I 



ION (00 
) ' l 1 



(29) 



Superscript K on various r matrices indicate that 
they act in Keldysh space. g$" A K {k) denote inverse 
retarded, advanced and Keldysh Green functions for 
non-interacting electrons in the graphene sheet while 
ga' A ' K {k) are the corresponding Green functions for lead 
a. For the graphene sheet, they are given by 

9o*(k) 



9${k) 
9o(k) 



u> — Ak + iS 
2iSK(uj). 



(30) 
(31) 



K 



Here, K{u>) = l+2n,i?(w) where tif(w) is the usual Fermi 
Dirac distribution function. S is an infinitesimal regu 
larization parameter. For the non- interacting case 
merely serves as a regularization for the Keldysh func- 
tional integral. Because a finite self-energy term is an- 
ticipated from the coupling of graphene electrons to the 
leads g$ can be safely omitted here (i.e. g^{k) w 0)^&. 



A. Integrating out the leads 

We now integrate out the leads degrees of freedom 
in order to obtain an effective theory only in terms of 
fields defined on the graphene sheet. The inverse re- 
tarded, advanced and Keldysh Green functions for the 
leads, 9a ' A ' K \ are those corresponding to free fermions, 
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and because the leads are always in thermal and chem- 
ical equilibrium the Keldysh Green function is strictly 
related to the retarded and advanced Green functions 
via the fluctuation-dissipation theorem (FDT). They are 
given by 



9a ( k ) = u ~ £ k,k s 
g%(k) = 2££tanh 



LJ__Pa 

2T 



(32) 
(33) 



Upon integrating over the leads, the resulting self-energy 
action becomes 



-X K (k)T?T«-X K (k)T»T?}lP k , (34) 



where 



£*(*) = £ 



dk z 



C 2 

Sq 



= -»r (35) 



£ A *(jfc), 



and 



E*(fc) 



a 

x<5(w - e k - e fc J 



2T 



= —2^ ^""^ r a tanh 



2T 



(36) 



Here, F Q = irpt^ measures the effective coupling strength 
between the layer and leads, and r = Tl + Tr. p is 
the lead density of states to tunnel into the layer as- 
sumed to be constant. The frequency- independent damp- 
ing coefficient, T, and the vanishing real energy shift 



that result from our assumptions indicate that the bath 
is treated as an Ohmic environmental. Combining the 
actions Eqs l26l34l we obtain the dressed inverse Green 
functions for electrons in the graphene sheet 



g R (k) = u- \ k + iT = g A *(k) 
g K (k) = 2^^^ Q tanl/ a; ~ /iQ 



2T 



(37) 
(38) 



The negative imaginary part of T, R (k) leads to an irre- 
versible damping in the time-dependent Green function 
G R (k, t). The damping term formally describes deco- 
herence suffered by a propagating electron wave due to 
incoherent escape and injection of electrons into and from 
the leads. 

At this point, it is convenient to shift the energy scale 
so that all energies are measured with respect to p = 
(Pl +Mj?)/2. This is equivalent to the following mapping 

OJ — > OJ — p 

Ak — * — p 

Pa -> V a /2, 

where Vl,r — ±V and V = pl — Pr- We assume V > 0. 
Following this choice the inverse retarded Green function, 
Eg 1371 remains invariant while Eg 1381 becomes 



(k) = 2i r Q tanh 



uj - V a /2 
2T 



(39) 



Using the dressed inverse Green functions defined in 
Egs l37l39l the effective action for the graphene sheet is 

S*f f = I - 2U [A* cl A g + A^A,/] , (40) 

Jk 

where the inverse Green function matrix 5^,7 1 is now 
given by 



/ g R (k) 


A 9 

-g R (-k) 


-5k 





g K (k) 










\ 


A* 





Hk 


Ki 













™5k 





g R (k) 


A, 








g K (k) 


A cl 







,9k 


a; 


-g R {-k) 








Ki 










A c; 








g A (k) 


A 9 
-g A (-k) 


"9k 







A* ; 


g K (k) 








A* 





9t 













A c i 


~5k 





g A {k) 


A 9 




V o 





A* ; 


g K (K) 





flk 


a; 


-9 A (- 


k)J 



(41) 



B. The Mean Field Equations 

In closed equilbrium, solutions to the mean-field 
gap and number equations on the honeycomb lattice 



have shown that while graphene exhibits a BCS-BEC 
crossover behaviour away from the Dirac point for in- 
creasing attractive interaction strength, u, superconduc- 
tivity in graphene at half-filling requires a finite attrac- 
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tive interactionJ^*i£. In this section, we derive the main 
results of our work which are the mean-field gap and 
number equations in the presence of leads and voltage. 
Solving these equations will allow us to study the effects 
of dissipation and nonequilibrium current on the gap as a 
function of attractive interaction strength u and filling n 
and compare the results to the equilibrium calculations. 
We begin by obtaining an effective theory for the s-wave 
order parameter alone by integrating out the graphene 
electrons. From Eq[3Dl we obtain 



dA* 



= 0. 



(43) 



A c ,=A : A o =0 



In principle, the action may be extremized with respect 
to A c i but the corresponding saddle-point will not be 
pursued here since it only gives a trivial relation. Eq]43l 
yields, 



iS£ / (A,A*) = Trln[- 



- 2iU (A* cl A q + c.c.) . (42) 



The gap equation 



° b eff 



OA* 



-iTr ■ 



A,=0,A e! 



A o =0,A c 



2A , s 

~u' (44) 



The self-consistent equation for the gap can be ob- 
tained from the following classical saddle-point equation 



This equation leads to the generalized nonequilibrium 
gap equation, 



J 



2A 



4Aa)^Q r « talm 



-V a /2 
2T 



[((w + £ k ) 2 + T 2 )((lu - £ k ) 2 + T 2 ) + 4A 2 | 5k | 2 ] 



[(^-£; + (k)) 2 + r 2 ][(c 



^-(k)) 2 + r 2 ][( w 
i 



£+(k)) 2 + r 2 ][(w + £_(k)) 2 +r 2 ]' 



(45) 



The spectra of the two bands are given by 



A 2 £±(k) = A k ±| 5k |, (46) 



and Ek = ^A 2 . + |gk| 2 + A 2 . After scaling all energies 
by bandwidth t and evaluating the cj-integral we obtain 



1 

u 



1 



2irN 



£[F„(S + (k))+^(E_(k))], 



(47) 



where 



F v (x) = - 
x 



tan 



tan 



and 



H±(k) = 



t 



U T a V 

t ' t ' t 



7 = 7l + 7_r denotes the sum of lead-graphene tunnel- 
ing rates scaled by t. Eqf?7] is the BCS gap equation 
in the presence of leads (7) and voltage (v) and is the 
nonequilibrium generalization of Eq.2 in refO Indeed 
when one takes the limit as 7 — > and v — > in Eql47l 
the equilibrium gap equation is recovered. 

At low-energies, excitations in graphene at or near 
half-filling are concentrated near two inequivalent Fermi 
points at the corners of the hexagonal Brillouin zone. In 
the vicinity of these points, we have 



Ai 



3t' 



fj, = m 



|.g± K +k| w v F \k\, (48) 



where vp — it/ 2 is the Fermi velocity and ±K = 
(±47r/3\/3, 0) are the locations of the inequivalent Fermi 
points. Within this approximation, the quasiparticle dis- 
persions, 3-t(k), become 



e±(k)«m±e S±(k) 



(49) 



where e = Uj?|k|. Noting that the area per lattice site 
is A/N = 3-\/3/4 the conversion from k-summation to 
e-integral is given by 



V 



3\/3 



N ^ 4tto| j 



ede, 



(50) 



The energy cut-off, set by conserving the total number 

of states in the Brillouin zone, is .33 in 

units of t. In the continuum limit, the gap equation then 
becomes 



1 



3^3 
8tt 2 v 2 f 



( D edc[f e (S + (k))+F (S_(k))]. 
Jo 



(51) 



2. The number density equation 

In equilibrium, the number density is computed us- 
ing a thermodynamic relation OTmf I 'dji — —N e . Out 
of equilibrium, the relation does not hold and the par- 
ticle density, n, must be extracted from one of the four 
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Kadanoff-Baym Green functions, G < usin g 22 i 23 

" = XE / G -A( fc )- ( 52 ) 

a labels the electron spin and A £ {A, B} labels the 
sublattice in which it propagates. In terms of Keldysh 
Green functions 1 ^, 




where G R ' A ' K (k) are the retarded, advanced and Keldysh 
Green functions for the graphene electrons. These 
Green functions can be obtained by inverting the matrix, 
<g- l (k), in EqHB We find that the form of the Green 
functions is independent of spin and sublattice, and the 
resulting number equation reads 



_ 47 r duj [1 - F{uj,v)]{cqlo & + c 5 w 5 + c 4 w 4 + c 3 w 3 + c 2 w 2 + C\LJ + c ) , 
H ~ * k J ^ + S +) 2 + + S -) 2 + 7 2 ][(" - 2+) 2 + 7 2 ][(^ - S-) 2 + 7 2 ] ' ( j 

F(ui, v) is the zero-temperature nonequilibrium electron distribution, and is given by 

An exact evaluation of the w-integral in Eq l54l is difficult. However, it can be done in the limit where the applied 
bias is assumed small compared to the bandwidth and the dampling coefficient, i.e. v <C min{l,7}. Computing the 
integral up to quadratic order in v the number density yields 

eae 



4t™ 2 Jo 

c (10 7 2 + S 2 ,) + ( 7 2 + S 2 h )( 7 2 + S 2 _)[2c 2 + c 4 (2 7 2 + E 2 + + S 2 _) + c 6 (10 7 4 + 6 7 2 S 2 _ + Bj + &y 2 E 2 + + 5 4 )] 

(7 2 + S 2 +)(7 2 + S 2_ )(1 g 7 4 + E 2 +(8l 2 + £ 2+ _ ^ + g3 (g^ + g2 & +)) 

2 



7T [(~\ - S 2 _) 3 + 8 7 2 S 2 ,(2 7 2 + ~l) - 8 7 2 S 2 _(2 7 2 + S 2 J] 

{tan~ 1 f" + ) 
E [ + 7> [c^l S 2 _ - 4 7 2 ) + c 3 (S 4 + 7 2 S 2 _ + 4 7 4 - ~ 2 + (~ 2 _ - 3 7 2 )) 

-c 5 [~l + 6 7 2 S 4 + + 9 7 4 H 2 + - S 2 _( 7 4 - 6 7 2 H 2 + + ~ 4 + ) - 4 7 6 ]] -(--+) 



In 



7 2 +S 2 



7 2 



|i) (2 Cl + e 3 (3 2 + + ~ 2 _ + 2 7 2 ) + 2c 5 (~ 2 - 2 _ - 7 2 ^ - 7 2 S 2 _ - 3 7 4 )) } 



+ {2X 1) ^( 7 2 + S 2 h ) 2 ( 7 2 + S 2 _) 2 ' ;+ 27r(7 2 + S 2 h ) 2 (7 2 + S 2 _) 2W ' 

(56) 

where x = 7l/ 7 and H± are given by EqJlSJ The coefficients Co, ... ,C6 are dependent on H±, £± and 7, and are 
defined as 

C6 = 1, 

"2 1 ~2 

c 4 = 3 7 - 



2 

c 3 = 2[e + ( 7 2 -S 2 _)+e-(7 2 -S 2 ,)], (57) 

c 2 = 3 7 4 + (El ~ 2 _) 2 + 7 2 (S 2 _ + S 2 + ) - ^ - =± 
Cl = e-(7 4 + 27 2 S 2 h +S 4 )+e + (7 4 + 27 2 S 2 _+S 4 J, 
co = i(S 2 _ +7 2 )(S 2 ,+7 2 )(5 2 + +S 2 _ + 27 2 ). 



I 

It can be easily verified that in the limit of 7 — > and (c.f. Eq.3 in refllq). The mean-field equations Eqs l51l56l 
v — > 0, Eq |56l rcduces to the equilibrium number equation 
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are the central results of this work. These equations will 
be analysed in the remainder of the paper. 



IV. RESULTS 

Our main focus will be on obtaining and analyzing 
gap phase diagrams in the parameter space of interaction 
strength (u) and number density (n) for various leads- 
graphene couplings (7l,7_r) and external biases (v). A 
previous work on closed equilibrium graphene^ revealed 
that, at half-filling, the superconducting instability of the 
semi-metallic phase requires a critical attractive interac- 
tion strength u c and, thus, the gap vanishes up to u c . 
Away from half-filling, the metallic phase is immediately 
unstable to superconductivity for arbitrarily weak attrac- 
tive interaction strength. As a result, the gap remains 
finite for any finite u and the system displays a typical 
BCS-BEC crossover behaviour. In this section we quan- 
titatively discuss the effects of dissipation and nonequi- 
librium current on the gap phase diagram by numerically 
solving the generalized mean- field equations, Eqs l51|56l 
The following sections will show that a dramatic mod- 
ification to the phase diagram is observed by the mere 
coupling of graphene to its environment, even in the ab- 
sence of nonequilibrium current. We find that the effects 
of external biases in addition to dissipation does not sub- 
stantially alter the qualitative features of the phase di- 
agram from the case in which the system is subject to 
dissipation alone. However, as the following sections will 
discuss, the application of an external bias leads to shifts 
in the metallic region surrounding half-filling which re- 
sult from voltage-induced changes in the graphene elec- 
tron density. The results presented here are applicable to 
the case of small biases (y <C min{l,7}); effects of large 
biases are not considered here. 



A. Finite lead-layer coupling 7 / 0, zero voltage 

(v = 0) 

First, we begin with the case in which the lead- 
graphene-lead heterostructure is in thermodynamic equi- 
librium. In particular, this is the situation where fJ-L — 
M-R = Mres, and in the long-time limit /jL sys = fi res is 
maintained. Here, electron tunneling processes between 
the central graphene system and the leads is providing a 
mechanism for decoherence for the particles in the sys- 
tem (7 7^ 0), but an external bias that explicitly breaks 
time-reversal symmetry of the heterostructure is absent 
(v = 0). Consider the case where the central graphene 
sheet is in a superconducting phase. Because of its cou- 
pling to the leads one can envisage a situation in which an 
electron that constitutes a Cooper pair escapes into the 
leads. Because the leads are assumed to be infinite the 
electron that has escaped the system is completely lost 
in the leads and as a consequence looses its coherence 
with its former partner. Although a different electron 
may enter the system from a lead within a time-scale of 



Ttun ~ 1/r, the electron will not necessarily pair with 
the widowed electron since it completely lacks coherence 
to do so. Because dissipation effectively acts as a pair- 
breaking mechanism we expect a suppression of the gap 
throughout the entire region of the phase diagram. 

FigEl plots the gap phase diagrams for various leads- 
graphene coupling strengths (7). FiglSJa) corresponds 
to the closed equilibrium case which has been obtained 
previously^. Figj5]^b) , (c) display the behaviour of the 
gap as 7 is increased. It is apparent from these plots 
that the suppressed region in the gap (dark blue region) 
grows as 7 is strengthened. Regions of large gap values 
corresponding to the region with large u also displays 
an overall suppression in the gap as 7 is increased. The 
qualitative features of the diagrams are consistent with 
the expectation described above. Let us now discuss the 
results more quantitatively. 



1. Half-filling (n = 1) 

For the closed equilibrium case at half-filling (7 = i> = 
and n = 1) the semimetal-superconductor transition is 
possible mainly because the divergent nature of the inte- 
gral on the right hand side of Eq[5T] is cured by particle- 
hole symmetry. When the integral is convergent, it is 
clear that a solution to the gap equation does not ex- 
ist for small u where becomes larger than the in- 
tegral. The value of the critical interaction parameter 
at which the transition occurs can be easily quantified. 
At half-filling the number equation, EqfSH is satisfied by 
m = 3t' — fi = 0, and thus at the critical point (n = 1, 
A = 0, and m — 0) the gap equation reads 



1 



3V3 



de 



1 



2.33 



(58) 



For any u < u c the equations cannot be solved with any 
real A and the system enters the semimetallic phase. In 
the presence of dissipation (7 > 0) the number equation 
is still solved by m = at half-filling, and the gap equa- 
tion at the critical point yields 



1 

Ur, 



3\/3 

Ml 



D 



detail 



3V3£> 
2^4 



e 

7 

7-D 



tan- 1 ( 7D 1 )-f ln(l+ 7D 2 ) 



,(59) 



where the reduced coupling strength is given by 7^ = 
•y/D. The integral on the right hand side of EqlSUl 
is convergent and, thus, tells us that the semimetal- 
superconductor transition exists in the presence of dis- 
sipation at half-filling. The behaviour of func- 
tion of jd is plotted in FigJU We see that the value of 
u c increases as 7 is increased. This is consistent with 
the above considerations from which we expect that a 
larger interaction parameter is necessary to achieve pair- 
ing since leads-induced decoherence generally suppresses 
superconductivity. The phenomenon can also be ob- 
served in Fig[5] where the apex of the blue region shifts 
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right for larger 7. The plots show that at 7 = 
verges to the closed equilibrium value of u c 
predicted by previous calculations. 



u c con- 
~ 2.33 as 




FIG. 4: The plot of critical coupling function of re- 

duced leads-graphene coupling 70 = y/D. 




1 2 3 4 5 

Interaction strength u [in units of t] 

FIG. 5: Plots of the BCS gap, A, in the parameter space 
of attractive interaction strength u and electron density n. 
The three diagrams correspond to different values of leads- 
graphene coupling strengths. In (a), the system is closed, i.e. 
7 = 0, while in (b) and (c) 7 — 0.1 and 7 — 0.2, respectively. 
As the coupling is increased, the blue region in the phase di- 
agram, where the gap is small, grows. In parts of the blue 
regions in (b) and (c) the gap is zero even for ti / 1, indi- 
cating that a metal-superconductor quantum phase transition 
emerges in the presence of dissipation. 



What is notable in Figj5]is the expansion of the blue re- 
gion, where the gap is small, as 7 is increased. The ques- 
tion is whether or not the typical BCS-BEC crossover 
behaviour observed in the closed equilibrium case is a 
correct physical picture away from half-filling for finite 7. 
The external baths acting as a pair-breaking mechanism 
makes the issue subtle. The pair-breaking perturbation 
in a superconductor with magnetic impurities has been 
shown 24 '- 5 to strongly suppress the transition tempera- 
ture of the superconductor. Therefore, when such pertur- 
bation is strong enough the gap may vanish completely 
and give rise to a metal-superconductor quantum phase 
transition at finite doping. The question of whether or 
not the gap vanishes away from half-filling depends on 
the convergence of the integral in the gap equation. At 
v = 0, the generalized gap equation becomes 



— oc / tat 
u Jo 



1 



— tan" 1 — +( 



7 



(61) 



We see that for any m (i.e. regardless of being at half- 




2 3 4 5 

Interaction strength u [in units of t] 

FIG. 6: The dark areas above show regions in the phase di- 
agram where the gap equation lacks a solution for any hnite 
A; the gap vanishes in these regions. As in Fig[5] the system 
is closed for plot (a) while 7 = 0.1 and 7 = 0.2 in plots (b) 
and (c), respectively. 



2. Away from half-filling (n 7^ 1) 

In the closed equilibrium case away from half-filling, 
m ^ and the critical point condition becomes 



1 



3V3 



cdc 



1 



1 



= 00. (60) 



The divergence of the integral results in a solution with 
A > for any small u > 0. This gives u c = implying 
that Cooper instability occurs for any finite u away from 
half-filling. Let us now investigate how this is modified 
when 7 is finite. 



filling or not), the integral is convergent because for any 
small 3±, which is the source of divergence, the arctan 
factor nullifies the divergence. This implies a finite u c 
at both half-filling and away from half-filling. Conse- 
quently, the system should undergo a superconductor- 
to-metal phase transition as the interaction parameter is 
lowered. Notice that the analysis above infers that the 
system will eventually enter the metallic phase as u is 
decreased for any density. 

Fig [6] explicitly shows regions in the gap phase diagram 
where the gap equation lacks a solution with any posi- 
tive A. The diagrams are plotted for the same values 
of 7 as in Figj5j The black regions are where the gap 
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equation is solutionless and represents a (semi)metallic 
phase. Clearly, as 7 is increased, the metallic region ex- 
pands. We find that the superconducting (white) and 
metallic (black) regions are separated by a second-order 
phase transition. 

The fact that the central graphene sheet becomes 
metallic away from half-filling certainly defies expecta- 
tions based on the intrinsic properties of graphene. For 
closed graphene the semimetal-superconductor quantum 
criticality emerges at half-filling because of the absence of 
electron-electron screening that results from the vanish- 
ing density of states at the Dirac point. As stated earlier 
in this section, the physics of the metallic region can be 
understood using the phenomenon of leads-induced de- 
coherence, and, thus, is not described by intrinsic prop- 
erties of the central graphene sheet. 



B. Effect of voltage, v ± 

So far, we have discussed the effect of leads-induced 
dissipation on the gap phase diagram in the absence 
of voltage. We now consider the effects of driving an 
out-of-plane charge current though the superconducting 
graphene sheet. Here, we are limited to the regime of 
small voltages, specifically v <C min{l,7}. As mentioned 
before, we assume v oc /il ~ M_r, >0 and allow the rela- 
tive strengths of the two couplings to the leads, 7l and 
7r, to vary. In the absence of current (v — 0), the gap 
equations depends only on the sum of these couplings 
7 = j L + "f R . But EqEH shows that in the presence 
of current (v 7^ 0) the number density now depends on 
these couplings independently and depending on the rela- 
tive strengths of these couplings the dominant correction 
term may change sign. The main qualitative modifica- 
tions to the gap phase diagram in the presence of finite 
voltage reflects the influence of this correction term. 

In the small voltage regime and for 7 < 1, the domi- 
nant correction term gives a correction of order 71* -C 1 
to the number density, which is of order unity. Because 
the modifications to the gap phase diagram due to volt- 
age is expected to be small the effect is more clearly seen 
by plotting the difference in gap values at finite and zero 
voltage. This is shown in Fig|7] The gap difference is 
plotted here for 2v = 7 = 0.2 in the vicinity of the 
apex region. In Fig[7Ja) Jl/jr — 4 while in FigJT^b) 
IlHr = 0.25. The plots reveal that the metallic region 
(blue region in FigJS]) shifts vertically away from half- 
filling. The figure shows that for Jl/jr — 4 the apex 
shifts up while for jl/jr = 0.25 it shifts down. Given 
M = (a*l + A*fl)/2 and v > the lowest order voltage 
correction in Eg 1561 tells us that the number density in- 
creases or decreases depending on the asymmetry of the 
lead couplings. If > Jr, n increases, while if jl < 1r, 
n decreases. The gap equation yields the largest value of 
u c given 7 and v when m — 0. Thus, the above observa- 
tion tells us that, for 7^ > Jr, m = is achieved not at 
half- filling as in the equilibrium case but at n > 1. This 
shifts the apex upward. The opposite occurs for jl < jr. 





(a) 








(b) 
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u/t 



FIG. 7: A cartoon plot showing the effect of voltage on the 
boundary of the metallic region. The dashed lines in both 
plots denote the boundary at v — 0. The shaded area is the 
metallic region after a steady-state bias is applied. In both 
plots, the applied boltage is v = 0.1. However, yt, > Jr in 
(a) while ~/l < 7h in (b). Essentially, the voltage-induced 
modification is to shift the metallic region to higher values in 
density or to lower values depending on the polarity of the 
voltage and the lead-coupling asymmetry. 

The nonequilibrium gap equation is convergent for all /x, 
thus, a metallic phase is once again expected at all den- 
sities. 



V. CONCLUSION 

In conclusion, we have theoretically studied the effects 
of dissipation and nonequilibrium drive on the properties 
of superconducting graphene. An external steady-state 
current was perpendicularly driven through the graphene 
sheet by attaching it to two leads which were equilibrated 
at two constant, but different, chemical potentials. The 
mean-field BCS theory of superconductivity on graphene 
was extended to the nonequilibrium situation by formu- 
lating the theory on the Keldysh contour. After obtain- 
ing nonequilibrium gap and number density equations 
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FIG. 8: A plot of ii c vs. v for a fixed fi. The plot line sepa- 
rates the metallic and superconducting phases of our system. 
Adjusting /i will tune the location of u° on the :r-axis but the 
general shape of the curve is not modified. 
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we studied the BCS gap as a function of attractive in- 
teraction strength u and electron density n for various 
lead-graphene coupling strengths 7 and voltages v. We 
have shown that dissipation results in a suppression of 
the BCS gap at both zero and finite voltages. We argued 
that the coupling of the graphene sheet to external baths 
acts as a pair-breaking mechanism because it causes an 
electron that constitutes a Cooper pair to escape into the 
leads. Once an electron leaves the scattering region, it 
looses coherence with its time-reversed partner and the 
destruction of the Cooper pair entails. 

A quantitative understanding of why the gap is signifi- 
cantly suppressed by dissipation can be gained by observ- 
ing how dissipation affects the gap equation. Recall that 
the BCS gap equation for an ordinary superconductor 27 
in closed equilibrium is given by 



this case, 



A = uTN(0) 



V" 2 n + A 2 



(62) 



iV(0) is the density of states at the Fermi energy, and 
u > is the attractive interaction strength. A general 
result for these ordinary superconductors is that the gap 
equation fEql6"2"]). and hence the gap, is unaffected by 
time-reversal-invariant perturbations. Take, for example, 
the influence of non-magnetic impurities on the supercon- 
ducting state. The gap equation obtained after invoking 
disorder-averaging and the Born approximation reads 



A = uTN(O) 



A 



(63) 



A 2 



where Co and A are frequency and order parameter renor- 
malized by the perturbatio n 26 i 28 i 29 , and N(0) is the den- 
sity of states in the presence of the perturbation. The es- 
sential point is that Co and A are related to their unrenor- 
malized counterparts by a common factor 77 — 77 (w n , A), 
i.e. 

Co = 1]LO, 

A = 77 A. 

Because this factor 77 cancels out in Eq[531 the gap equa- 
tion remains invariant and leads to the result that the 
gap is unaffected by non-magnetic impurities 30 . 

Imagine now that a pure ordinary superconductor is 
coupled to an external bath in equilibrium. The Nambu- 
Gorkov equations can be straightforwardly derived for 



(ico n + isgn(io n )T - £ k ) G + AF = 1 (64) 
(iio n + isgn(uo n )T + f k ) F + AG = 0. (65) 

where the ordinary and anomalous Green functions are 
given by 

G(k,co n ) = - f dr(r TCk , T (r)c kiT (0))e iw " T 



dr (r r c k , T (r)c_ ka (0))e^ r . 



We immediately see from Eqs l64l651 that 10 and A scale 
asymmetrically, namely, 

r 

Co = rjLO A = A; 77 = 1 + 1 r . (66) 

\v n \ 

Here, T is the rate at which electrons decay into the 
bath. The asymmetry in the renormalization of 10 and 
A (Eq l66|) greatly affects the gap equation, Eq[63j and 
shows how dissipation can affect the gap significantly. 
This is consistent with the qualitative argument given 
above. 

The emergence of the metal-superconductor quantum 
phase transition in the graphene subsystem at both zero 
and finite voltages gives rise to the possibility of induc- 
ing the phase transition using external bias. While fixing 
the average chemical potential fj, to some value, v can 
be changed by adjusting /_t£ and fin symmetrically about 
fi. u c is obtained from the gap equation in this situation 
by fixing A = and fi to some value. FigJS] shows a 
generic plot of u c as a function of voltage. If the inter- 
action strength,?/, of the system is at u — u*, then for 
v < v(u*) the system will be metallic. However, when 
v is increased and passes v = v(u*), the system will be- 
come superconducting. vP c can be tuned by adjusting 
the average chemical potential fi. It is clear from Eg 1561 
that when the average chemical potential fi is fixed, the 
electron density can change as a function of voltage. 
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